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ABSTRACT 

We develop an approximate analytical solution for the transfer of line-averaged radi- 
ation in the hydrogen recombination lines for the ionized cavity and molecular shell 
of a spherically symmetric planetary nebula. The scattering problem is treated as a 
perturbation, using a mean intensity derived from a scattering-free solution. The an- 
alytical function was fitted to Ha and H/3 data from the planetary nebula NGC6537. 
The position of the maximum in the intensity profile produced consistent values for 
the radius of the cavity as a fraction of the radius of the dusty nebula: 0.21 for Ha and 
0.20 for H/3. Recovered optical depths were broadly consistent with observed optical 
extinction in the nebula, but the range of fit parameters in this case is evidence for a 
clumpy distribution of dust. 

Key words: stars: AGB and post-AGB - stars: evolution - Infrared: stars — (ISM:) 
planetary nebulae: general - ISM: jets and outflows - ISM: molecules - 



1 INTRODUCTION 

Shells with internal cavities are found in almost all PNe 
(Balick & Frank 2002). Such shells and cavities are related 
to the shaping of PNe by stellar winds. Low and intermedi- 
ate mass stars lose their mass at a very high rate, forming 
a circumstellar envelope. Fast stellar winds from the cen- 
tral stars of PNe overtake the slower, denser, AGB ejecta, 
resulting in interaction of the fast and slow stellar wind com- 
ponents (Kwok, Purton & Fitzgerald 1978). The fast wind 
sweeps the slow wind, and leaves a cavity inside. Cavities 
are also found in supernova remnants, (Dwek et al. 1987; 
Lagage et al. 1996) but these structures are more compli- 
cated (Bilikova et al. 2007). Shell-cavity dimorphism is also 
found in LBV stars and in AGN. 

In one of the well studied bipolar planetary nebulae, 
NGC 6537, the core consists of bright arcs tracing a shell sur- 
rounding an elongated cavity (Matsuura et al. 2005). Arcs 
are bright in Ha and other recombination lines, and extinc- 
tion maps derived from Ha and H/3 suggests the presence 
of dust grains in the arcs. Matsuura et al. (2005) suggested 
that little dust exists in the cavity. Deriving accurate dust 
density from Ha and H/3 line maps is not straightforward, as 
the light from the central star is scattered within the arcs. 
In these arcs, both dust grains and gas are mixed together. 
Therefore, we developed a radiative transfer code to resolve 
cavity and shells for PNe, which includes scattered light in 
a shell. 



A self-consistent radiative transfer code for this mixture 
has been developed by Ercolano, Barlow & Storey (2005). 
They have used a Monte Carlo method. Here, we use an 
analytical solution, with some approximations, and concen- 
trate on the configuration of cavities and shells around them. 
The aim of the paper is therefore to obtain detailed physi- 
cal insight through a simplified analytical model, which can 
complement the more complex information from numerical 
solutions. 



2 THE MODEL OF NGC6537 

In this section, we describe the physical and radiative trans- 
fer models that are used in our analysis of NGC 6573. Al- 
though certain parameters are definitely specific to this ob- 
ject, taken from Matsuura et al. (2005), much of the model 
is generally applicable to any planetary nebula of similar 
geometry, and in a similar state of evolution. 



2.1 Physical Model of the Nebula 

We assume that the nebula NGC6537 is accurately spheri- 
cally symmetric. The layout of the various radial zones of the 
object is summarised in Fig. 1. The central star is a white 
dwarf of negligible solid angle, both from an observer's point 
of view, and from the modeller's: no rays in the radiative 
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4" 

1 .5 Extinction 

Figure 1. A diagram of the planetary nebula model. The spher- 
ically symmetric nebula is centered on the star, A, of negligible 
solid angle. This is surrounded by the cavity, B, and molecular 
shell, C. An infinitessimal volume, D, of the cavity is marked at 
radius, r. Optical radiation, E, escapes from the shell after suf- 
fering some degree of extinction. The particular values shown (in 
magnitudes) apply to NGC6537. 



transfer model are considered to start or end on the stellar 
surface. 

The central star emits sufficient vacuum ultraviolet ra- 
diation to ionize a surrounding cavity. Ions in this cav- 
ity undergo frequent recombination and photoionization cy- 
cles, and it is the main source of Ha and H/3 line radi- 
ation. For NGC6537 specifically, the observed flux ratio 
Fhp/Fhoi = 2.79, where these fluxes are averages over 
the appropriate spectral-line bandwidth. From this ratio, 
we adopt radially constant values of the electron temper- 
ature, T e — 1.5 x 10 4 K, and electron number density, 
n e = TO 4 cm" 3 . 

Surrounding the Hn cavity is a neutral shell, composed 
mainly of molecular material. This shell contains dust, which 
we assume, for the case of NGC6537, to have a mass fraction 
of 0.01. Most of the shell material is, of course, in the form of 
molecular hydrogen. The shell produces approximately two 
magnitudes of visual extinction in NGC6537, derived from 
the measured excess, E(Ha — H/3). The shell is spatially 
resolved, and the extinction varies from 1.5 magnitudes at 
the centre to 2.2 magnitudes at 4" off centre. 

At present, we leave the form of the density profile as 
an unknown function in both the cavity and shell zones of 
the nebula. 



2.2 Radiative Transfer Model 

The volume element D in Fig. 1 has volume SV. If we let this 
be unit volume, then the radiated power per unit frequency 
and solid angle in the Ha line is 

j v {Ha) = ^^-A 3 2n 3 (r)(t> v ,32, (1) 

47T 

where ^32 is the 'laboratory' line-centre frequency of the Ha 
line, and A32 is the Einstein coefficient for spontaneous emis- 
sion. The subscripts refer to the principal quantum num- 
bers of atomic hydrogen, so nz(r) is the number density 
of H-atoms in the upper state of the transition. With an 
isothermal approximation, the lineshape function, <f> v ,32 is 
not a function of radius, but may be a function of direction 
if there are significant radial motions in the cavity or shell. 
Given these definitions, a similar expression may be written 
down for the emission coefficient of H/3: 

j v (Hp) = l ^lA A2 n 4 {r)(p vA 2. (2) 

We assume complete velocity redistribution, so that absorp- 
tion lineshape functions are identical to those specified above 
for emission. 

The radiative transfer equation is 

h-VI v = -[n c (r) +a c (r) + K„(r)]I v 

+ K c {r)B„(T{r))+j„(r) + <j c (r)J v {r), (3) 

where ft is a unit vector along the ray of specific intensity 
I v . Opacity is provided by radius-dependent absorption co- 
efficients, n v , k c , for the line and continuum, respectively, 
and the continuum scattering coefficient, <j c . Kirchoff 's Law 
allows the continuum emission from dust to be written as 
K c (r)B„(T(r)), whilst the line emission coefficient is j v (r), 
and will be one of the specific forms in eq.(l) or eq.(2) for the 
appropriate transition. The final term on the right-hand side 
of eq.(3), equal to a c (r)J v , is the scattering integral assum- 
ing isotropic, elastic scattering, and J„ is the angle-averaged 
intensity. 

We expand the left-hand side of eq.(3) in spherical polar 
coordinates (for example Peraiah (2002)). This introduces 
H = cos# where 9 is the angle between the direction of ray 
propagation and the radial direction. On the right, we com- 
bine the continuum absorption and scattering into an extinc- 
tion coefficient, X c ( r )i an< A assume that the line absorption 
is negligible when compared to that in the continuum. The 
transfer equation with these modifications is 

+ a c (r)J„+j„(r) (4) 

We now integrate eq.(4) over the spectral- line bandwidth ap- 
propriate to either the Ha or the H/3 line. We assume that 
this bandwidth is adequate to cover all the line radiation 
in the two hydrogen lines studied, regardless of all Doppler 
shifts within the source. In this connexion, we note that the 
velocity extent of the line due to the bulk motion of expan- 
sion is of order 20kms _1 , whilst the thermal Doppler width 
is 21.4T 4 1/ ' 2 kms -1 , with T4 equal to the kinetic temperature 
in the ionized cavity in units of 10 4 K. An unknown micro- 
turbulent width must be added in quadrature to the latter 
figure. It is therefore reasonable to suppose that even the 



© 0000 RAS, MNRAS 000, 000-000 



Radiation Transfer in PNe 3 



extreme red- and blue-shifted portions of the line are signifi- 
cantly blended by the thermal and microturbulent lineshape. 
In other words, the combination of a low terminal expansion 
velocity and a large thermal plus microturbulent line width 
allows us to neglect velocity field induced Dopplcr shifts that 
would otherwise complicate the analysis considerably. Let 
the spectral-line bandwidth be Av, and the line-integrated 
intensity is given by 



/A1//2 
-An/2 



I v dv 



(5) 



A similar equation to eq.(5) relates the angle-averaged in- 
tensities J and J„. If we assume also that the functions x°, 
k c , B V (T) and a c vary only very slightly over Av, these 
functions can be removed from the filter integral when it is 
applied to eq.(4). The result is 



dl 

^ dr 



(1-H 2 ) dl 
r cfyi 



+ 



-X c (r)I + Avn c {r)B v {T{r)) 
a c (r)J + j(r), (6) 



where j(r) = J^^ 2 jvdv. The only frequency-dependent 
part of j„ is the appropriate lineshape function (see eq.(l) 
and eq.(2)). We assume that the filter width is sufficient for 
the normalisation condition of the lineshape to hold, so that 
the integral j(r) in eq.(6) is either, 

j(Ha) = ^A 32 n 3 (r) (7) 

or the equivalent expression for H/3. 

Finally, we re- write eq.(6) in separate forms appropriate 
to the cavity, and to the shell, respectively. In the cavity, we 
make the approximation that there is no dust, and that other 
processes, such as free-free and bound-free emission, make 
a negligible contribution to the radiation flux within the 
filter bandwidths used, for example Matsuura et al. (2005). 
The cavity therefore acts effectively as a pure source of line 
radiation, and the radiation transfer equation in this zone 
is, 

dl (1 - n 2 ) dl 

+ - " ' 

or r 



dfi 



= j(r). 



(8) 



By contrast, the shell, which is rich in dust, and where al- 
most all the hydrogen is molecular, does not emit any Ha or 
H/3 line radiation, is a site of continuum absorbtion, scatter- 
ing and thermal emission by dust. The shell transfer equa- 
tion is therefore, 



dl (1 - n 2 ) dl 



+ AuK c {r)B v {T{r)) + o c {r)J (9) 



We now proceed to solve eq.(8) and eq.(9) in turn, coupling 
the solutions via boundary conditions. The solution is pre- 
sented in several stages, and we outline these briefly here. 
In Section 3, we solve eq.(8) for the specific intensity of line 
radiation in the cavity, and this solution is angle-averaged in 
Section 4 yielding a mean intensity in the cavity as a function 
of radius. This quantity is not required for later calculations, 
so Section 4 is only included for completeness. In Section 5, 
we solve eq.(9) for radiative transfer in the shell: initially the 
specific intensity is obtained for dust extinction only, but a 
scattering source is then re-introduced as a perturbation. 




Figure 2. A ray enters the cavity, of radius R, at the top of the 
figure along path s in direction n. The ray is initially at an angle 
ui to the radial unit vector, r. This relationship is modified along 
s as f changes direction: at an arbitrary point at radius r the ray 
unit vector ft and the radial unit vector cross at angle 6. Note 
that when s is smaller than the distance to the mid-point of the 
chord, 8 becomes an obtuse angle. 



The cavity and shell solutions, for the specific intensity of a 
ray that may pass through both zones, are combined in Sec- 
tion 6. Section 7 is devoted to developing a mean intensity 
in the shell, as a function of radius, by angle-averaging the 
unperturbed solution from Section 6. This angle-averaged 
intensity is used to develop explicit forms for the scattering 
perturbation in Section 8. The remaining sections are con- 
cerned with observational properties of the solution, and fits 
to HST data for NGC6537. 



3 THE CAVITY SOLUTION 

The layout of a typical ray passing through the cavity zone 
of the nebula is depicted in Fig. 2. The model cavity has no 
inner radius. All rays enter the cavity from the surrounding 
shell with some boundary value of the specific intensity. This 
value may differ from ray to ray, and, at present, we leave 
this value as an unknown. The emission coefficient depends 
on the radius only through the number density of H atoms in 
the upper level of the appropriate Balmer line from eq.(7). 
We make the approximation that these number densities are 
constants in the cavity, on the grounds that the cavity gas 
is in pressure equilibrium because of its high sound speed, 
and that we have already assumed a constant temperature 
in this region. We ignore the effect of overpressured regions 
near the edge of the cavity where the ionized medium is 
juxtaposed against the surrounding shell (Perinotto et al. 
2004). Therefore, we can re-cast eq.(8) as 
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dl (i-^)dl 

A 4 n T 



dr 



dfi 



= Jo, 



(10) 



where jo is the constant emission coefficient. The left-hand 
side of eq.(10) is the spherical polar expansion of the dot 
product, h • VJ, but an alternative representation is as the 
derivative, dl/ds, taken along the ray characteristic. An al- 
ternative way of writing eq.(10) is therefore, 



dJ _ . 
ds ^° 



which has the simple linear solution, 
I{s) = I + j s, 



(11) 



(12) 



where To is the unknown intensity with which the ray enters 
the cavity from the shell. The problem now consists only 
of re-expressing the distance along the characteristic, s, in 
terms of the independent variables r and /i = cos 9. To ob- 
tain such a relation, we apply the sine rule to the triangle 
in Fig. 2, which is bounded by r, s and R, yielding, 

r - R - 8 (13) 
sin a; sinfl sin(a; + 9) 

Now both R and u are constants, so it follows from the first 
equality in eq.(13) that 



r sin 9 = R sin u — K, 



(14) 



is a constant along the characteristic. Applying the cosine 
rule to the same triangle yields a quadratic equation in s. 
The geometry of Fig. 2 requires the positive root, so that 
the distance along the characteristic is, 

Nl/2 



s = r cos 9 + (R z 



Vsin 2 9f 



(15) 



Using the result in eq.(14), and the definition of fi, we obtain 
s = r\i + (R 2 — K 2 ) 1 / 2 . The constant expression under the 
square root reduces to J? cos a;, so the characteristic distance 
is s — rfi + R cos ui. The intensity along the characteristic is 
therefore, 



I(r,fi) = I + jo(r/i + -Rcoscj), 



(16) 



which is the required cavity solution for constant emission 
coefficient. It is trivial to prove, from eq.(16) that it satisfies 
eq.(10). 



4 MEAN INTENSITY IN THE CAVITY 

To obtain the angle- averaged intensity, J{r), we need to in- 
tegrate over a number of different rays, which all meet at 
one point, where the radius, r, is a constant. At this point, 
9 is the final angle along its path, and can be used in the 
required solid angle integral. However, now uj is a variable, 
as we are integrating over many rays, and u needs to be ex- 
pressed in terms of 9, or /i. The situation is summarised in 
Fig. 3. Results from the cavity solution, which still apply to 
Fig. 3, are that sin a; = (r/R) sin 9 and s = r cos 9 + R cos w. 
We use the first of these relations to eliminate w from the 
second. The expression for s can in turn be used to write 
the cavity solution, eq.(16) in the form, 



I(r,n) = I (n)+j (rn+[(R 2 



2 2il/2 



+ r n ] 



(17) 



noting that n - 
side of eq.(17). 



cos 9 is the only variable on the right-hand 




Figure 3. The angle averaged intensity, J, is to be computed at 
point P, at radius r from the centre. The average is over rays 
(examples marked with arrows) which make angles, 6, with the 
radius vector at P, where 6 may be in the range to ir. A partic- 
ular ray is shown entering the cavity at angle ui, in direction n, 
and passing through a distance s in order to reach P. 



The solid angle integral for the mean intensity, with the 
integration over the azimuthal angle already carried out is 
f i 



J(r) = 



J(r, n)dn. 



(18) 



After substitution of eq.(17) into eq.(18), and a partial eval- 
uation of the resulting integral, we find, 



J(r) 



Jo + -^jor 



(a + /J 2 ) 1 'dn 



(19) 



where a = [(R/r) 2 — 1], and Jo is the angle average of I . 
The result is a standard integral in terms of the arcsinh 
function, and after converting this to logarithmic form, the 
angle-averaged intensity is 



J(r) - Jo + If 



r + R 



VR 2 



(20) 



We note that the apparent infinity in eq.(20) at r — disap- 
pears when this equation is replaced by a suitable expansion 
for the case r <C R- The small radius form is, 



J(r) -> Jo + (l/2)ioi?(2 - r 2 /R 2 



(21) 



We plot, in Fig. 4, the function Q = (J(r) — J )/(j R) from 
eq.(20), as a function of the dimcnsionless radius x = r/R. 



5 THE SHELL SOLUTION 

We re-write the shell transfer equation as 



-X c (r)I + f(r), 



(22) 



di | (1-mV J. 

dr r <9/x 

where f(r) is an arbitrary function of the radius, which in- 
corporates the continuum emission and line scattering. We 
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5.1 Shell Solution with Power-Law Density 




0.2 0.4 0.6 0.1 

Dimensionless Radius 



Figure 4. A dimensionless form of the angle- averaged intensity 
in the cavity, [J(r) — Jo]/ (jgR), see eq.(20), plotted as a function 
of dimensionless radius, r/R. Equation 21 is used for the point at 
r/R = 0. 



change the left-hand side of eq.(22), as in the cavity solu- 
tion, to describe a solution along the ray characteristic. The 
equation can then be put in standard first-order linear form 

as 



dl 
ds 



+ X C (r)I = f(r(s)). 



(23) 



Equation 23 may be integrated by the standard method of 
integrating factors. The boundary condition at s — is that 
the specific intensity enters the shell from interstellar space 
with the intensity Ibg, assumed equal to a typical value for 
Galactic starlight. The extinction coefficient at this same 
position is zero, regardless of its variation within the shell. 
With this condition imposed, eq.(23) has the formal solu- 
tion, 

1(a) = I BG e- -lo xC(r(s ' ))ds ' 

+ J /(»•(«')) e*pj-jf X c (r(a))da|d S '. (24) 

The ray distance, s, is related to the radius r and direction 
cosine fi by a modified form of the relation which appears 
in eq.(16). If the shell has an outer radius R2, then 



s — rfi + R2 coscj, 



(25) 



where, for the moment, we ignore the presence of the cavity. 
With this same caveat, r sin 8 = R2 sin a; is a constant along 
a given ray if uj is the angle with which the ray enters the 
shell from interstellar space. Defining K2 = R2 sinoj, we can 
express jj, in terms of the radius as \i = 

(l~K'i/r 2 ) 1/2 , which 
may be substituted into eq.(25) to yield a relation between 
r and s with no other variables involved. This relation, with 
the radius as the subject, is 



r = (s — 2si?2 cosw + R 2 ) 



2x1/2 



(26) 



The expression for the radius in eq.(26) can be used to ex- 
pand the radius in terms of s in eq.(24), provided that func- 
tional forms for x c ( r ) an d f(r) are known. 



To progress beyond eq.(24) we require some analytical ap- 
proximation for the behaviour of the radius-dependent func- 
tions. We start with the extinction, \ c ■ We assume that this 
depends on the radius only through the number density of 
the shell dust, so that the functional form is the same for 
both the absorption and scattering contributions. Further, 
we assume a density dependence which follows the inverse 
square behaviour of the singular isothermal sphere. There- 
fore, we suppose that the extinction in the shell behaves as 



X c (r) = x C {R)(r/Ry 



(27) 



where X C {R) ls the maximum value of the extinction coeffi- 
cient, found just outside the cavity boundary. With the help 
of eq.(26), we can write the extinction as a function of s 
rather than r. Integrals of the type which appear in eq.(24) 
can now be written in the form, 



/ 



X c (r(s))ds = R 2 X C (R) 



ds 



(28) 



Introducing the new variables x = s — R2 cos u> and a = 
R2 sin a;, eq.(28) may be written as the standard integral, 



J x C (r(s))ds = R\<(R) J 

which has the solution (as an indefinite integral), 



X c (r(s))ds = R x ( R ) arctan 
R2 sin uj 



R2 sin u> 



(29) 



(30) 



Equation 24 contains two definite forms of eq.(30). The first 
of these has limits of to s, and can be re-written, with 
the help of addition formulae for the arctangent from Grad- 
shteyn & Ryzhik (1965), as 



X c (r(s'))ds' 







{tv} + arctan ^ 



R2 



,(31) 



where the n is to be included only for ray distances such 
that s > R2/COSU1, and the group /3 = R 2 x c {R)/R2- The 
second definite integral form from eq.(24) has the limits s' 
to s, and it is convenient to write it without combining the 
arctangents, since the part with s can be removed from the 
integral over s' . Overall, the shell solution can now be writ- 
ten, 



I(s) = Ibg exp 



{it} + arctan 



(R2/ s) — COS UJ 



+ e si 
where 

u(s) = arctan 



^P (s ')exp(^ ( ,,<. 



S — R2 COS UJ 



(32) 



(33) 



R2 sin uj 

and the functional form of / remains to be determined. 
5.2 Scattering as a perturbation 

If we ignore continuum emission in the shell, the unknown 
function, / in eq.(32) reduces to line scattering alone, which 
depends on the angle- averaged mean intensity, J{r). It is 
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then possible, in principle, to average eq.(32) over solid an- 
gle, leading to an integro-differential equation for J(r). How- 
ever, owing to the difficulty in attempting to solve such an 
equation analytically, we resort to the simpler procedure of 
treating the scattering term as a perturbation. For this ap- 
proximation to be very good, the scattering contribution 
to the extinction should be small, so that c c (r)/x c ( r ) is a 
small parameter for all radii. This is unlikely to be true for 
real dust. For example, silicate dust modelled by Ossenkopf, 
Henning & Mathis ( 1992) has an optical efficiency for scat- 
tering which is consistently ~2-3 times that for absorption 
over optical wavelengths. However, we still adopt the pertur- 
bative approach as the only viable method of obtaining an 
approximate analytical solution. In this procedure, we first 
solve the radiative transfer equation in the shell, asuming 
that line extinction is the only contribution to the right- 
hand side of eq.(22). This means that we can set / = in 
eq.(32), and take as the zero-order solution in the shell, 



I(r, n)= Ibg exp 



[to 



+arctan 



ril-fi 2 ) 1 / 2 

H ^ 2 ) 1/2 [^+(^-r 2 (l-M 2 )) 1/2 ' 
r(l-/^)- M (fl2-r2(l-^))i 



/2 



,(34) 



which is just the first term of eq.(32), with 9 and u> elimi- 
nated in favour of /i = cos 6, and with the help of the re- 
lation sin lo — (r/ife) sin # (see eq.(14), which still holds in 
the shell). It is tedious, but straightforward, to show that 
eq.(34) is indeed a solution of eq.(22), for the case where 
f(r) = 0. Equation 34 cannot be used alone, as along any 
path through the nebula, a ray may encounter a sequence of 
shell, cavity, and again shell conditions. We therefore need to 
determine some functional form for Ibg for the case where 
a ray leaves the cavity and re-enters the shell, forming a 
combined solution along a ray. 

As scattering is now to be treated as a perturbation, 
it can be computed from a mean intensity which is derived 
from the zeroth-order combined solution. An equation for a 
perturbation in the intensity, SI, can be constructed from 
cq.(9) by expanding the specific intensity as I — J» + SI, 
and then subtracting off the equation in the zeroth-order 
estimate, The result is, 



d(SI) (1 - fj, 2 ) d(SI) 



- X c (r)SI + * c J(r), 



(35) 



dr r d/i 

where we have assumed that continuum emission makes a 
negligible contribution, so that f(r) = a c J(r). So, mathe- 
matically, the equation in the perturbation may be treated 
in the same way as the shell equation, eq.(22). The formal 
solution for 81 therefore looks like eq.(24) with the scattering 
term replacing the unknown function /. The perturbation 
solution is therefore given by 



61(e) 



= SI(so)e 



-h(s) 



I a c (r(s'))J(r(s'))e h(a ' ) 

J Sn 



ds , 



(36) 



where the function h is given by eq.(30). It is important to 
note that there are three possible versions of eq.(36): A ray 
which avoids the cavity has <5/(so) = and a lower limit of 
so = on the remaining integral. A ray which enters the 
cavity has two shell segments. In the first, 8 1 (so) = and 
so = 0, but the upper limit is the value of the path where 



R 



/or 




Figure 5. The complete ray path is shown divided into three 
segments: between P and Q a distance s p is covered in the shell, 
having entry angle u>. The ray enters the cavity at Q, with angle 
u>', and covers the distance s' between Q and Q' before exiting 
into the shell. The point M is the mid-point of s', and the closest 
approach of the ray to the central star (at O). Once in the shell 
again, the ray reaches the general point X, at radius r, where 
its path makes an angle, 9, with the radial direction. The ray 
direction vector and radial direction vector (at various positions) 
are marked with bold and light arrows respectively. In progressing 
from Q' to X, the ray covers the distance s s h. e ii- The ray finally 
leaves the nebula at point Z, where r = fu- 



tile ray enters the cavity. The perturbation is then constant 
in the cavity, forming a finite value of /(so) for the second 
shell segment, where s is now the path length where the 
ray leaves the cavity. 



6 THE COMBINED SOLUTION ALONG A 
RAY 

Given the approximations made in the previous section, it 
is now possible to combine the cavity and shell solutions for 
a single ray. The geometry of the combined solution is set 
out in Fig. 5. As we are ignoring, as this stage, rays which 
pass solely through the shell, each effective ray has three 
segments: the first is an absorptive passage through the shell. 
However, if the input background at r = R2 is very weak, 
we can ignore this segment, and treat the specific intensity 
of the ray at the cavity boundary, r — R as zero. 

The second segment of the ray path, between points 
Q and Q' in Fig. 5, passes through the cavity. Here the 
specific intensity is assumed to increase through spontaneous 
emission in the line, according to the cavity solution, eq.(16). 
In particular, we want an expression for the specific intensity 
at Q', where the cavity solution becomes the input value for 
the shell solution, which takes over as the third segment of 
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the ray path, which continues until point Z is reached, and 
the ray passes into the vacuum. 

At the point X in Fig. 5, which is as radius r, the ray 
has travelled a distance s s hell through the shell, and it makes 
an angle 9 with the radial vector. At this point, the solution 
will be the shell solution, with an input specific intensity 
from the cavity. 

Analysis of the triangle MOQ shows that the cavity 
segment of the ray path has length s' = 2Rcosuj' . The total 
distance from point P to point M is, from the triangle MOP, 
equal to sm = R2 cos oj. Therefore, the 'pre-cavity' distance 
(from P to Q) is s p — R^cosu — Rcosuj'. The additional 
distance along the ray from M to X is s M i — r cos 9 = r/i. 
The total distance along the ray from P to X is therefore 
s — R2 cosuj + r/j,, just as in the shell solution, eq.(25). From 
these results, we find that the distance s s heii in Fig. 5 is given 
by 



Sshell =rfJb— R COS Uj', 



(37) 



and by applying the sine rule to the triangle OPQ, the angles 
uj' and u> are related by 

sin a/ = (R2/R) sinw, (38) 

which can be used to eliminate u' from eq.(37), yielding 

Sshell = ryi — [R 2 — R% sin 2 cj] 1/2 . (39) 

The limiting (maximum) value of uj, for a ray which just 
enters the cavity, occurs when ui' — 7r/2, when sinuwi = 
R/R-2- 



6.1 Input Solution to the Shell 

As we are assuming that the specific intensity at point Q 
in Fig. 5 is effectively zero, the input solution for the shell 
is just the cavity solution evaluated at Q'. When this is 
done, the cavity solution, eq.(16) applies, with uj 1 replacing 
uj. We take this solution with Jo = 0, r = R and /1 = cos a/, 
obtaining I(R,uj') = 2joRcosuj' . With the help of eq.(38), 
this expression becomes 



I(R,uj) = 2j [R 2 



,1/2 



(40) 



which is the background for the solution in the shell beyond 
point Q'. 



6.2 Combined Solution in the Shell 

As the limits of integration are modified for the combined 
solution, we generate the shell solution from the indefinite 
integral, eq.(30). The upper limit is the distance along the 
ray path, s, but the lower limit is now, from Fig. 5, s' + s p — 
R2 cosw + [R 2 — R 2 sin 2 uj] 1 ^ 2 . The result for the extinction 
integral is 



/: 

J s' - 



X c {x)dx = 



R 2 x c (R) 

R2 sin uj 



arctan 



S — R2 cos u 



R2 sin uj 



arctan 



(R 2 



R\ s^ujf' 2 



R2 sin uj 



(41) 



When eq.(25) has been applied, and eq.(41) has been re- 
expressed entirely in terms of fi, the result is, 



X c (x)dx = 



I { 

(1-^)1/2 \ 



arctan 



arctan 



(1 - M 2 )!/2 



(l-(r/i?) 2 (l~M 2 )) 1/2 



(42) 



(r/flXl-M 2 ) 

The first arctangent can be converted to a simpler form, as 
an arcsine, via the relation 



arctan 



(1 -x 2 ) 1 / 2 



(43) 



which is taken from Gradshteyn & Ryzhik (1965). The sec- 
ond arctangent in eq.(42) can be written as the reciprocal 
of the form which appears in eq.(43) by defining the new 
variable q = r(l — fi 2 )/R. An additional relation linking arc- 
tangents (Gradshteyn & Ryzhik 1965) for the case where 
x > 0, which is true of q, is 



arctan(l/a;) = 7r/2 — arctan a;, 



(44) 



and this relation allows us to use eq.(43) directly on the sec- 
ond arctangent in eq.(42). A third relation from Gradshteyn 
& Ryzhik, arcsina; + arccosa; = it/ 2 allows us to re- write 
eq.(42) as, 



X c (x)dx = — 



13 



,2W2 



[arcsin q — 9] . 



(45) 



It is now a simple matter to write out the combined solution 
for a ray by inserting eq.(45) into the shell solution, eq.(24), 
with Ibg given by eq.(40). With /3 and q fully expanded in 
terms of 9, the result is 

(R 2 /R) 2 sin 2 uj] 1 ' 2 



I(r,9) 



2joR[l 



exp 



-Rx c (R) 

(r/R) sin 9 



R 



,(46) 



which is valid for R < r < i?2- For plotting, we re- write 
eq.(46) in the form 



I(x,a) = 2j R(l 



2\l/2 

■ a ) exp 



t . .a 
- — (arcsin a — arcsin — 1 
a x 



,(47) 



an equation in two parameters, a = (R2/R) sino;, and the 
optical depth parameter, r = Rx c {R)- The latter parameter 
becomes equal to the radial optical depth of the shell in the 
limit where R2 R- The independent variable is x — r/R. 
This last definition implicitly introduces a third parameter, 
the upper bound on x, equal to R2/R. We note that eq.(47) 
reduces to a sensible limiting form, that is 



I(x,a = 0) = 2j -Rexp[-r(l - R/r)], 



(48) 



for the case where a = 0. In Fig. 6 we show the function 
U — I(x,a)/(2joR) for five values of a between its mini- 
mum value of zero, and maximum of 1. We set the optical 
depth parameter to be r = 1.38 (see Section 10) to agree 
approximately with the visual extinction of 1.5 magnitudes 
near the centre of NGC6537, and we plot the dimensionless 
radius out to x = 3. 



7 MEAN INTENSITY IN THE SHELL 

The approximation discussed in Section 6.1 - that rays that 
do not cross the cavity contribute zero specific intensity - 
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1.5 2.0 2.5 

Dimensionless Radius 



3.0 



Figure 6. The function U = I(x,a)/(2joR) (see eq.(47) plotted 
as a function of dimensionless radius x = r/R for values of a = 0.0 
(solid line), 0.4 (dotted line), 0.6 (dashed line), 0.8 (simple chain), 
0.99 (complex chain). 



excludes all rays with negative values of cos#. Rays with 
non-zero specific intensity are limited to a subset of those 
with positive values of cos#, more precisely to those with 
sin 9 < R/r. With the additional restriction that any radius 
where J(r) is computed has R < r < R2, the situation here 
is similar to that shown in Fig. 3. The function to be aver- 
aged is the combined solution specific intensity in the zeroth- 
order (no scattering) approximation, given by eq.(46). The 
above considerations allow us to write a formal integral for 
the mean intensity, recalling that r sin 9 = R 2 sincj: 

/>arcsin(i?/r) 

J(r)=j R [1- (r/R) 2 sin 2 9] 1/2 



x exp 



Rx c (R) 

(r/R) sin 9 



/ r sin 9 \ 
\~R~) ~ 



sin 6>d<9 (49) 



We simplify eq.(49) via a series of substitutions. The 
first of these is to let the impact parameter be p — 
(r/R) sin 9; we also define the new constant parameters, 
7 = R\ C (R), and p = R/r. These definitions transform 
eq.(49) to, 



J(r) 



1-P 2 
1 - P V 



1/2 



x exp ^ [arcsinp — arcsin (pp)] f dp, (50) 



noting that p and p are always < 1. Given this condition, we 
can expand the inverse sines in eq.(50) in terms of Gauss hy- 
pergeometric functions (for example, Abramowitz & Stegun 
(1965)). For complex argument z, 



arcsin;? = zF(l/2, 1/2; 3/2, z 2 ), 



(51) 



where the power series forming the Gauss hypergeometric 
function, F, is absolutely convergent within, and on, the unit 
circle for the arguments in eq.(51) (Gradshteyn & Ryzhik 
1965). The substitution of eq.(51) into eq.(50) has the ben- 
eficial consequence of cancelling p within the exponential, 
and leaving p 2 everywhere except for the product pdp. This 



suggests the substitution x = p , yielding, 

2 /■! / -, \ 1/2 



J(P) 



joRp 







1-x 
1 - p 2 x 



X exp {- 7 [F Q, i; §, x) -p F [\, I; | p\)]} dx(52) 

At this point, we expand the hypergeometric series which 
appear in eq.(52). For the case here, where the first two 
arguments are the same, the power series (Abramowitz & 
Stegun 1965) is 



F(a,a;g,z) 



a 2 (a + l) 2 z 2 
+ gV. + 5(5 + 1)2! 
a 2 (a + l) 2 (a + 2) 2 z 3 



g(g+l)(g + 2)3\ 



+ ■ 



(53) 



in the case of general a, g, and for the specific case of a = 1/2 
and g = 3/2, we find, 



F(l/2,l/2;3/2, Z ) = l + 



3z 2 



5z A 



40 + 112 + 



(54) 



Substitution of eq.(54) in eq.(52) leads to the expression 

2 /■! / -, \ 1/2 



J(P) = 



joRp 



X 



1 — p 2 x 



exp{- 7 [(l - p) 



+ (l-p 3 )x | 3(l-p^)x 2 | 5(1 -pV | 



6 



40 



112 



,(55) 



which is as far as it is possible to proceed without making 
further approximations. 



7.1 Approximate Forms 

Although the series in eq.(55) is convergent for all relevant 
values of x = [(r/R) sin#] 2 , convergence is rather slow when 
x is close to 1, the condition for a ray in glancing contact 
with the cavity. However, the contribution of large x to the 
integral is limited by the term in front of the exponential 
which tends to zero as x — > 1 (for the particular case of 
p — 1, see below). The physical reason for this is that rays 
at large angle have relatively short paths through the cavity, 
and correspondingly small specific intensities on entry to the 
shell. 

Although values of p which are close to 1 (positions close 
to the shell/cavity boundary) can force the leading term in 
eq.(55) to 1, this is not a problem because, in this situation, 
the argument of the exponential approaches zero, the whole 
exponential to 1, and integration of the leading term alone 
is sufficient for moderate accuracy. We evaluate eq.(55) to 
two levels of accuracy in x. In the zero-order approximation, 
we abandon all but the first term in the series, but this does 
not depend on x, and can therefore be removed from the 
integral, leaving, 



J(p) 



joRp 2 



-7(1"P) 



a* 



p 2 X 



1/2 



dx. 



(56) 



On evaluation of the integral in eq.(56), and reverting to the 
original radius variable, r = R/p, we obtain the zero-order 
approximation to the mean intensity: 



ji \ joR - 
J(r) = —e 



7(l-fl/r) 



R 2 



2rR 



(57) 
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Figure 7. The function Q = J(x)/(joR) plotted as a function 
of dimensionless radius x = 1/p = r/R using the zcroth-order 
approximation, eq.(57), (solid line), the first order approximation, 
eq.(59), (dashed line), and exact integrals from eq.(50) (diamond 
symbols). 



We note that at the cavity/shell boundary, where r = R, the 
reduction of eq.(57) agrees with the cavity solution, eq.(20), 
evaluated at the same radius. Both equations yield J{R) = 
j R/2. 

It is also possible to obtain an analytic integral for a first 
order approximation, keeping the term in x in the series in 
eq.(55). The integral in the first-order case is a standard 
form in Gradshteyn & Ryzhik (1965). With the parameters 
specific to the current problem, the integral is 



l-x 
1 — p' 2 x 



1/2 

exp 



- 7 (1 - p 3 )x 



dx = B(l, ■ 



x*i(i,|,|;p 2 .«) 



(58) 



where z = —7(1 — p 3 )/6. The first function, B, on the right- 
hand side of eq.(58) is an Euler beta-function, whilst the 
second is a confluent hypergeometric series in two variables, 
p 2 and z. Further standard formulae allow the beta-function 
B(l, 3/2) to be expressed as a ratio of factorials that reduces 
to 2/3, so that the mean intensity to first order in x is 



j(p) = JoV c -7d-P) 



1 I I 2 

' 2' 2' P '' 



-7(1 -P 3 ) 



(59) 



We plot, in Fig. 7, the mean intensity in the shell as 
a function of dimensionless radius, 1/p, computed from the 
exact integral, eq.(50), and both approximations (eq.(57) 
and eq.(59)). For the value of the optical depth parameter 
used (7 = 1.38, see Section 6.2), we see that even the zeroth- 
order formula is an excellent approximation. The graphs in 
Fig. 7 should be seen as an extension of that in Fig. 4 to 
values of r/R > 1 and to Q < 0.5 in the case where the 
input background is negligible. 



8 PERTURBATION SOLUTION IN THE 
SHELL 

We now proceed to solve eq.(36) with the mean intensity in 
the scattering term given by eq.(57): it is too complicated 



to use the more accurate eq.(59). There are various forms of 
cq.(36) which should be used for the appropriate zone of the 
source. For a ray which penetrates the cavity, we initially 
consider the case where this ray has entered the shell, but 
has not yet reached the cavity. Assuming a negligible input 
intensity from the vacuum, we have the boundary condition 
that so = and 51(0) = 0, and therefore the problem for 
this ray, and zone, reduces to solving the integral, 



(60) 

which appears 
(61) 



81{s) = e- h(s) ( o- c (r(s'))J(r(s'))e h(s,) ds'. 
Jo 

The most useful form for the function h(s') 
in the integrating factor, for this zone is 

, = fly(fl)[7r-a,-fl( 8 ')] 
R2 sin to 

which is derived from eq.(31), with a lower limit of and 
an upper limit of s' , followed by transformations similar to 
those used in working from eq.(43) to eq.(46). Note in par- 
ticular that for a ray in this case, uj is an acute angle, but 
9 is obtuse, and the respective limits of these angles for a 
radial ray are zero and 7r. 

We substitute for the scattering coefficient, a c (r(s')), by 
assuming that it has the same 1/r 2 functional dependence 
as the absorption (see Section 5.1). The mean intensity is 
given by eq.(57). It is perfectly possible to express all these 
functions in terms of the distance, s', along the ray, but the 
integral in eq.(60) appears easier when working in terms of 
radius. Letting the integral be we have 

a c (R)R 2 jo exp [e(ir - w - a)] 



dp exp [— e(arcsin ap — op)] 



(1 



z 2 p 2)l/2 



(1 



2p 



1 + P 



(62) 



where p = R/r as before, a — (R2/R) sin u, e = j/a, and 
the limits are given by 

Po = R/R2 (63) 
and 

P = R/[s 2 -2sR cos cj + Rl} 1 ' 2 , (64) 

which reduces to 1 for a ray entering the shell from the 
vacuum, and reaching the edge of the cavity, where r = R. 
If we divide eq.(62) by the top line, which is independent of 
p, we can carry out an integration by parts. The integrated 
part is 



_ f exp [— earcsin(ap)] 

V -J (l_ a 2 p2)1 /2 



exp [— e arcsin(ap)] 



(65) 



which removes the problematic square-root term. The result 
of this integration is, 



-e(a — ap+arcsin(ap)) 



(1 



2p 



In 



1 + P 



f 

JR/R 2 



g(p) exp { — e (a — ap + arcsin(ap))} 1 



R/R 2 
(66) 



where the function g(p), now entirely composed of loga- 
rithms and rational functions of p, is defined by 



g(p) = 2ea ■ 



2i + 
P 



ea 1 
1 1 5- + tap 

p p 2 



In 



1 + P 



(67) 
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and Vl/j is given by 



4ea* 



a c (R)R 2 j e<T 



(68) 



Note that the argument of the exponential in eq.(68) is dif- 
ferent from that in eq.(62) because the original integral, that 
is the lower line of eq.(62), was multiplied by a factor of 
2eae~ ea in order to obtain ^ j in eq.(66). 

To calculate an approximation to the integral in eq.(66), 
we expand the argument of the exponential, noting that 
ap < 1. The first term in the power-series expansion of the 
arcsine cancels with op, and the next term is in (ap) 3 which 
we ignore, such that 



-e(a + arcsin(ap) — ap) 



-ea(l + (ap) 3 /(6a)) 



(69) 



which is independent of p, and can be moved outside the in- 
tegral, which has now been reduced to the integration of the 
function g(p) from eq.(67). The integration of g(p) breaks 
down into six integrals, five of which can be solved in terms 
of elementary functions and the sixth in terms of the Rie- 
mann $(z,s,v) function (Gradshteyn & Ryzhik 1965). A 
full form for the solution, appears in Appendix A. 



8.1 Rays which avoid the cavity 

A ray which does not enter the cavity obeys the same prop- 
agation equation, eq.(60), but has a different upper limit on 
the distance, s. Instead of integrating to the cavity bound- 
ary, we integrate to the mid-point, where the radial vector 
in perpendicular to the direction of the ray. At this point, 
the value of s — ifecosaj, and the corresponding radius 
is 7?2sinai. We can therefore still use eq.(62), eq.(63) and 
eq.(64), but noting that the expression in eq.(64) now re- 
duces to P — R/(B,2 sin w), rather than 1. The integrations 
and approximations which follow still apply, providing we 
use the new value of P. 



8.2 Outward bound rays 

Provided that we are interested only in emergent rays, we 
can use the symmetry of the nebula to simplify matters 
greatly here. The integral along the ray, V&j, is indentical, 
whether we are integrating inward to the midpoint (or cav- 
ity) , or outward again towards the edge of the nebula. These 
integrals become formally the same because when we con- 
vert from s' to r as the integration variable, we must use 
different roots of the expression 



' o i/2 n 2 • 2 sl/2 

s = H.2 cos u) ± (r —Pi sm lo) ' , 



(70) 



with the negative root required for the inward segment, and 
the positive root holding for the outward ray. The result 
is that the integrals over radius are identical. The form of 
eq.(60) which applies to outward rays is 

8I(s) = 8I(s ) + e- h(s) I a c (r{s'))J(r{s'))e h< - a ">ds', (71) 

where so is now either the midpoint, for a ray which avoids 
the cavity, or the outward cavity boundary. We assume in 
eq.(71) that there is negligible scattering within the cavity. 
The input, 51 (so), is just the inbound solution, evaluated at 
the midpoint or the inbound cavity boundary. If we define, 




1.0 
Ray Angle 



Figure 8. The function U = SI(s)/(2joR) for complete path 
lengths s = s(i?2) plotted as a function of entry angle into the 
nebula, ui. The optical depth parameter is \ C (R)R = 1.38 as 
before, and the ratio of the shell to cavity radii is R2/R = 3.0. 
The ratio of the scattering coefficient to the extinction coefficient 
is a c {R)/x c {R) = 0.25. 



C = a c (R)R 2 j /(4ea), 

then we can re- write eq.(71) as 

§J(s) — | e e arcsin [(-H2/i-(so)) sinu] 



+ 



arcsin[(H2 / r ( s )) sin w 



} 



(72) 



(73) 



noting that any constants in the exponentials have been can- 
celled between the e~ h term outside the integrals, and the e h 
inside. The equivalence of the integrals, tyj , has been used, 
from the symmetry argument above. We can now write down 
two versions of eq.(73). The first is for a ray which avoids the 
cavity, where r(s ) = i? 2 sinoj and r(s) = R 2 on emerging 
from the cavity. In this case, 

5I(s) = C,^ 3 e elT,2 (l + e -e(-/ 2 -«0). (74) 

For the other case, where the ray crosses the cavity, the first 
value changes to r(so) = R; the second remains the same, 
so we have, 



5I(s) = C*,(e £ 



(75) 



where the integral fyj is given in Appendix A (with appro- 
priate limits). 

In Fig. 8, we plot the emergent specific intensity for 
scattering, U — 8I(s)/(2joR) over the full range of entry 
angles to the nebula. The value of s is the largest path length 
through the shell material for each angle. The same optical 
depth parameter for extinction, 1.38, has been used as in 
Fig. 6 and Fig. 7. The ratio of the outer (shell) radius to the 
inner (cavity) radius is 3.0: this has not been used as a formal 
parameter before, but the same value was also adopted for 
the above figures. The ratio of the scattering coefficient to 
the extinction coefficient has, rather arbitrarily, been set to 
0.25. For the perturbative method to be accurate, this value 
should really be small, but is unlikely to be so for typical 
dust models (see discussion in Section 5.2). This parameter 
anyway acts as a simple scaling factor which does not change 
the shape of the function in Fig. 8. The normalising factor 
of 2joR is the same as in the other specific intensity plot, 
Fig. 6. 
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The marked cusp in the curve in Fig. 8 is at the ex- 
pected angle, arcsin(i?/i?2), where the ray switches from 
paths which enter the cavity to paths which do not. The 
path at this angle also corresponds to the maximum dis- 
tance travelled through the shell material. The solution at 
smaller angles than that corresponding to the cusp comes 
from eq.(75); at larger angles, eq.(74) has been used. 



9 SOLUTION FOR EMERGENT RAYS 

To use the symmetry property of the integral tyj, we have 
already evaluated the perturbation, 81 as an emergent quan- 
tity above. In general, the complete solution, for a given 
emergent ray, is given by, 

Iem{s) = I(s)+SI( S ), (76) 

where s = 2i?2coso;. For a ray which crosses the cavity, 
I(s) comes from a form of eq.(46) where r = R2 and 6 — u, 
whilst SI(s) is given by eq.(75). Overall, 

J o-' D/1 „2\l/2,— efarcsin a — lj) . /-,t, / e arcsin a . ^£U}\ I *7*7\ 

lem = 2jon(l-a ) e '+QVj{e +e ).(77) 

A ray which does not cross the cavity has only the 51 con- 
tribution, and is given by eq.(74) without modification. 



9.1 The observer's view of the model 

To an observer with perfect angular resolution, an emergent 
ray of given exit angle, u>, characterises a circular strip of a 
spherical surface. This strip has area 27ri?2 sin it, where u is 
the polar angle measured from the line of sight towards the 
limb of the nebula. The observer cannot see this surface as 
a whole, but only a 2-D projection of it. The projected area 
of the strip is 

dA± = 2-kR% sin u cos udu (78) 

It is straightforward to see that for any given strip, the po- 
lar angle u is identical to the entry /exit angle of the ray, w, 
which has allowed values between and 7r/2. The observer 
with perfect angular resolution will therefore be able to pick 
out a small piece of a given (projected) strip, and will mea- 
sure the associated specific intensity, I em , as a brightness. 
For a nebula of radius R2, at a distance d from the observer, 
the measured brightness at an angle from the centre of 
the nebula, as seen by the observer, corresponds to exit an- 
gle lo = arcsin(0d/i?2). The brightness itself can be found 
by subsitituting this angle into either eq.(77) or eq.(74). 

The other extreme observer's view is that of a telescope 
with a beam that is much larger than the nebula. In this 
case, the specific intensity cannot be measured directly, and 
a flux, averaged over the whole object, is obtained instead. 
To a very good approximation, this flux is given by 

2-kr 2 n^ 2 

F — ^ 2 2 /(w)sinwcoswda;, (79) 

where, as previously, the functional form of is taken 

from eq.(77) or cq.(74), depending on whether or not the 
ray at angle to traverses the cavity. 



10 FITS TO NGC6537 

The function summarised in eq.(76) and eq.(77) was fitted 
to observational data describing the brightness variation of 
the Ha and H/3 spectral lines as functions of angular posi- 
tion across the nebula. The fits were made with respect to 
four variable parameters in the theoretically-derived func- 
tion: x — R/R2, the ratio of the cavity radius to the overall 
radius of the nebula; S = o c j\ c , the ratio of the scattering 
to extinction coefficients in the continuum; r = Rx c (R), the 
optical depth parameter, and an overall intensity-axis scale 
factor, Y, allowing a fit to normalized observational data. 
The optimum fit for each observational data set was taken 
to be that with the minimum value of the \ -statistic, 

2 1 ( h - Ie m (smLUi,X,S,T,Y)\ 2 

x -N=Z2^[ a, J' (80) 

»=i v 7 

for a data set with N entries of the form (sinoj^, /;), and the 
four free parameters introduced above. Details of the treat- 
ment of the observational data, and of the fitting function, 
are described below. 



10.1 Data from NGC6537 

Data were extracted from Hubble Space Telescope Ha and 
H/3 images (Matsuura et al. 2005). Seven straight-line slices 
were taken through the nebula, each at a different angle on 
the sky. For each slice, data were recorded for both Ha and 
H/3, and for each slice and spectral line, data were organ- 
ised into pairs consisting of a coordinate position along the 
slice, and a corresponding specific intensity. The following 
operations were applied to convert these data into a form 
suitable for fitting: For each slice and line, two pixel posi- 
tions were found for the emission peaks, corresponding to 
the intersections of the slice with the cavity boundary. The 
coordinate origin was then shifted to the mid-point of the 
peak positions, and the absolute value of the slice coordinate 
taken, transforming the data to a pair of radial brightness 
profiles measured from an origin at the approximate centre 
of the nebula. An approximate correction for the aspherical 
nature of the nebula was imposed by re-scaling the radial 
axis such that the origin to peak distance was the same for 
all radial profiles. The outer radius of the nebula was taken 
to be the smallest maximum value on the new scale, and the 
remaining profiles truncated to the same distance. Finally, 
the new radial coordinate was re-scaled once more to the 
range 0.0 — 1.0, with 1.0 corresponding to the now common 
outer radius. We note that the H/3 data covered a larger an- 
gular range on average than the Ha data, leading to a larger 
value of R2 in H/3 by a factor of 1.51, and a different radial 
scaling for the two lines. 

For individual radial profiles, the specific intensities 
were scaled such that the peak brightness was set to 1.0. 
Mean profiles for each spectral line were also constructed by 
averaging over all the un-normalized individual brightness 
profiles (on the normalised radial scale for each line defined 
above), and then normalising the averaged brightnesses to 
a peak height of 1.0. 

In the averaged profiles, the standard errors resulting 
from the averaging process far exceeded the formal mea- 
surement errors in the observational data. The former were 
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calculated at the peak position to be 0.098 for Ha and 0.054 
for H/3, as fractions of the mean peak height. The latter 
were only of order 1.7 x 10~ 4 and 3.0 x 10 -4 , respectively. 
For the normalised data sets, with peak heights of 1.0 we as- 
sumed absolute uncertainties at all radial points to be equal 
to those at the respective peaks: 0.098 for Ha and 0.054 for 
H/3. The same absolute errors were assumed for the fits to 
the individual radial profiles. Values of N — 4 were 75 for 
Ha and 120 for H/3. 



10.2 Parameter ranges 

The parameter x has the formal range 0.0 — 1.0, but a much 
smaller realistic range can be selected for any radial pro- 
file by observing the position of the scattering peak (see 
also Fig. 8). Both the optical depth parameter, r, and the 
scale-factor, Y, have the possible range 0.0 — oo. The most 
problematic parameter is S, the ratio of the continuum scat- 
tering to extinction coefficients. The scattering has been 
treated as a perturbation in Section 8, so strictly speak- 
ing only the range 0.0 — 0.1 is available to S. However, from 
dust models, we expect the ratio of the scattering to ab- 
sorption cross-sections to be of order 2 — 3 (corresponding 
to S = 0.666 — 0.75) in the optical region, as already dis- 
cussed in Section 5.2. We have therefore fitted each profile 
twice: once for a true perturbation, with S limited to the 
range 0.0 — 0.1, and once with the range 0.0 — 0.666, with 
the upper limit dictated by computed dust parameters. For 
each spectral line, we have fitted the averaged profile, and 
one selected individual profile. The fitting parameters for 
Ha and H/3 have been taken to be entirely independent in 
this preliminary study. 



10.3 Results of fitting 

In Fig. 9 we show the results of the fits where the allowed 
range of S is 0.0 — 0.666. The upper graphs are fits to Ha 
data, and the lower graphs, to H/3. For each line, the left- 
hand panels are fits to the averaged profile, whilst the right- 
hand panels show fits to a selected individual profile. Note 
that the individual H/3-fit (bottom right) is the only case in 
which the best fit fell in the perturbation range. A general 
difficulty with fits to averaged data is the broadening of the 
scattering peak in the averaging process, a consequence of 
the asphericity of the nebula that has only been approxi- 
mately mitigated by the rescaling operations described in 
Section 10.1. A full table of the fitting parameters and qual- 
ity estimates appears in Table 1. 

Figure 10 shows the fits in which S was restricted to the 
perturbation range, 0.0 — 0.1. Only in the case of the fit to 
the selected H/3 profile was the best overall fit found in the 
perturbation range, and this fit has already been shown in 
Figure 9. We note that this particular fit has rather peculiar 
parameters (see Table 1), with a much larger optical depth 
parameter (and much smaller scattering parameter, S) than 
the others. 

The easiest parameter to compare with observations is 
x. Values of x = R/R2 are markedly different for the two 
spectral lines, but this is simply a consequence of the greater 
angular extent of the data (larger R2) for H/3. When cor- 
rected to the value of R2 for Ha, by multiplying by 1.51, 



Table 1. Values of the fitting parameters, x,S,T and Y (see Sec- 
tion 10), that gave the best value of x 2 , as defined by cq.(80), 
for the associated data set. Type 'Avg' refers to averaged data, 
type 'Ind', to an individual profile. The symbol S' has the value 
3.2 X 10~ 4 . Note that values of x for H/3 should be multiplied by 
1.51 to place them on the same radial scale as those of Ha. 



Line/Type 


Smax 


X 


S 


T 


Y 


x 2 


Ha/Avg 


0.666 


0.2326 


0.666 


1.665 


0.372 


0.491 


Ha/ Avg 


0.100 


0.2120 


0.100 


2.413 


0.099 


1.423 


Ha/Ind 


0.666 


0.2060 


0.585 


2.850 


0.160 


0.384 


Ha/Ind 


0.100 


0.2060 


0.100 


3.175 


0.785 


0.434 


H/3/Avg 


0.666 


0.1320 


0.300 


1.602 


0.750 


0.755 


H/3/Avg 


0.100 


0.1325 


0.100 


1.995 


1.818 


1.079 


H/3/Ind 


0.666 


0.1346 


S' 


6.193 


26.10 


0.146 



the values of x for H/3, in the order they appear in Table 1, 
are 0.1993,0.2001 and 0.2032. These values are then con- 
sistent with those found for Ha. For the purposes of the 
optical-depth comparison (see below), it is useful to calcu- 
late uncorrected values of 1 — x for both lines, yielding 0.79 
for Ha and 0.87 for H/3. 

We now consider the results for the optical depth pa- 
rameter. In Section 2.1, the observed magnitude drops across 
the nebula, owing to extinction, ranged from 1.5 to 2.2 mag- 
nitudes, corresponding to intensity ratios of 0.251 and 0.132 
along a line of sight. The respective optical depths are 1.38 
and 2.03. We note that the fitted values of r in Table 1 are 
of the optical depth parameter introduced at the end of Sec- 
tion 6.2, and that this is related to the optical depth of a 
radial ray passing through the full thickness of the nebula by 
r r — t(1 — x) (see eq.(48) with the limiting value of r — > %)• 
Fitted values of r from the averaged data of 1.67 for Ha and 
1.60 for H/3 with the full range of S then yield respective op- 
tical depths of 1.32 and 1.39. If the perturbation restriction 
is enforced, the optical depths are 1.90 and 1.74, results that 
are reasonably consistent with the observational values Note 
that considerably different optical depths are recovered from 
the selected individual profile, suggesting a locally clumpy 
medium. 



11 CONCLUSION 

An approximate analytical function has been derived that 
predicts the brightness of a spectral line, within some filter 
bandwidth, as a function of angular distance from the cen- 
tre of a cavity/shell planetary nebula. The scattering part 
of this function is derived as a perturbation on the radiative 
transfer solution. This function has been fitted to observa- 
tional data from NGC6537 in the Ha and H/3 lines. The 
best-fit optical depth parameter in the model gives reason- 
able agreement with the observed optical extinction in the 
nebula. However, best-fit values for the scattering parameter 
are mostly too high to be consistent with a true perturba- 
tion. Consistent values were also found for the sizes of the 
cavity in both lines when the outer radius for Ha was ap- 
plied to both spectral lines. The mean of the values of x 
from Table 1 was 0.21, and the mean of the adjusted values 
for H/3 was 0.20. 

The analytical function may also have applications to 
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sin CO sin CO 

Figure 9. Best fit profiles (solid lines) following the functional forms in cq.(76) and eq.(77) to observational data from NGC6537 (dashed 
lines). The left-hand panel uses the observational data averaged over 14 radial profiles, whilst the right-hand panel uses just one individual 
profile. For each panel, the upper graph is for the Ha line, and the lower, for H/3. 



other sources with a hot, largely ionized, interior, sur- 
rounded by a shell rich in dust. Examples include post- 
AGB stars, supernova remnants, symbiotic stars, quasars, 
and highly-magnetized planets in addition to the planetary 
nebulae discussed here. 
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APPENDIX A: THE INTEGRAL OF G(P) 

The solution shown here is for any ray. For a ray which 
enters the cavity, the upper limit of these integrals should 
be set to P = 1. For a ray which avoids the cavity, the upper 
limit should be R/(R2 sinoj). The integral of g(p) in eq.(67) 
breaks down into six integrals. The first two, trivially, are 

h = 2ea ( dp = 2ea(P-R/R 2 ) (Al) 
J R/R 2 

and 

h= [ dp/p=2\n{R 2 P/R) (A2) 

The third can be integrated by parts, with the logarithm as 
the differentiated part. The result is, 




= (1 + P)ln(l + P) + (1-P)ln(l-P) 
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The fourth integral is the most problematic, as the re- 
sult cannot be expressed in terms of elementary functions. 
If the integral is split, so that 



h 
ea 



Jr/r 2 P V-PJ 
p 



- s 



1-P 



dp 



-f 

J R/R 2 



ln(l - p) 



dp, (A4) 



then both the expressions on the bottom line of eq.(A4) 
conform to a standard integral (Gradshteyn & Ryzhik 1965), 
which is their no. 2.728 (version 2). The solution is written in 
terms of the function $(«, s, v), which is discussed in detail 
in Section 9.55 of Gradshteyn & Ryzhik (1965). The final 
result is 



ea{P[$( 
R_ 

R 2 



-P,2,1) + *(P,2,1)] 
*^' 2 ' 1 > + *^' 2 ' 1 )]} 



(A5) 



The fifth integral can be solved in a similar way to the 
simpler 73, integrating by parts with the logarithm as the 
differentiated part. The integrated part of the result can be 
solved by an expansion in partial fractions, and the overall 
expression is 



h = 



1 



In 



R/R 2 



1+P 



dp 



R \R 2 -RJ \ 



Rl 



R 



R? 



+ 21nP - (1 + 1) ln(l + P) + (I - 1) ln(l - P^A6) 

The final integral can be solved in a broadly similar manner 
to I 5 , and the solution is 

.h, | } dp 

<R/R 2 



h 
ea 



2{ 1 -Rl) ln {R^R) +P -R- 2 
P 



+ 



1 



(1 



(A7) 



2 V ' Vl + P> 

The integral ^j, in its approximate analytical form, is 
now given by eq.(66), with its lower line replaced by the 
expression, 

6 

e- Ea ^/ fc , (A8) 
fe=i 

where the Ik are the integrals in eq.(Al) to eq.(A7). The 
necessary approximation to the exponential is discussed in 
Section 8 (see eq.(69)). 
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